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Asymmetric patchy particle models have recently been shown to describe the crystallization of 
small globular proteins with near quantitative accuracy. Here, we investigate how asymmetry in 
patch geometry and bond energy impact the phase diagram and nucleation dynamics of this family 
of models. We find the role of the geometry asymmetry to be weak, but the energy asymmetry 
to markedly interfere with the crystallization thermodynamics and kinetics. These results rational- 
ize the success and occasional failure of George and Wilson's proposal for protein crystallization 
conditions as well as provide physical guidance for developing more effective protein crystallization 
strategies. 

PACS numbers: 87.15.ak, 82.70.Dd, 87.15.nt 



Protein structures from X-ray and neutron crystallog- 
raphy are key to studying elementary biological mech- 
anisms [ll |2|. The difficulty of obtaining diffraction- 
quality crystals thus severely limits our understanding 
of living systems. From a physical point of view, iden- 
tifying successful crystallization conditions is equivalent 
to determining a protein solution phase diagram. Crys- 
tallizing a protein indeed typically involves placing a 
drop of protein solution near a high-salt aqueous buffer 
that drives the vapor diffusion of the solvent away from 
the drop. The non- volatile solutes then steadily concen- 
trate and, if the initial conditions are properly chosen, a 
protein crystal forms [T. The limited usefulness of ex- 
isting physical descriptions and of knowledge-based ap- 
proaches ^, however, leave a vast space of solution 
conditions to be experimentally screened. A material 
understanding of protein assembly is thus essential to 
developing more effective crystallization strategies. 

Soft matter descriptions of protein assembly based on 
particles with isotropic, short-range attractive interac- 
tions [5]-[7] - as suggested by early structural studies 0- 
[TO] - provide some conceptual guidance. They identify 
the region between the solubility line, above which the so- 
lution is stable, and the liquid-liquid critical point, below 
which the system gels [11 , as the "crystallization gap" 
where crystal assembly is possible. This simplistic pic- 
ture is, however, unable to reproduce many experimental 
trends [12-16 . The introduction of bond directionality 
in symmetric "patchy" models aims to better represent 
the effective protein-protein interactions that drive their 
crystallization [171 IH] • Yet the most commonly studied 
versions of these models have symmetric and interchange- 
able patches, which are atypical of real proteins [H [191 [20] 
and insufficient to describe the assembly of even the sim- 
plest of globular proteins [21, 22 . In this letter, we thus 
investigate the role of patch geometry and bond energy 
asymmetry on the phase diagram and assembly dynam- 
ics of a family of protein models. Note that this ad- 
ditional anisotropy 'direction' also complements earlier 
experimental proposals for introducing richness to soft 
matter assembly [23 and opens the door to assembling 



more complex structures in systems such as DNA-coated 
colloidal particles f24H26]. 

Model description - We consider patchy models hy- 
brid between those of Sear [27 and Kern-Frenkel [28] . 
Hard spheres, whose diameter a sets the unit of length, 
are decorated by three pairs i of patches that only in- 
teract with each other. Each pair's properties are in- 
dependently tuned [29 . The short radial extent of the 
square- well attraction, = l.la [30 , and its surface 
coverage measured by the semi-opening angle of its con- 
ical segment, 5i = cos~^(0.89), are chosen to be typical 
of protein-protein interactions [20l [22] By contrast, the 
patch orientation and bond energy are randomly cho- 
sen under the constraint that the lattice formed by simply 
bonding the patches is that of the most commonly ob- 
served in monomeric protein crystals, the orthorhombic 
P2i2i2i [311 E2]- This lattice's three non-intersecting 
two-fold screw axes guarantee a high number of rigid- 
body degrees of freedom with minimal symmetry con- 
straints. We characterize the patches with energy and 
geometry asymmetry parameters 

respectively, where li are the eigenvalues of the inertia 
tensor of the collection of equal- weight patches |29l |33] . 
Both C and 7 G [0, 1], where corresponds to an equal 
energy distribution (61=62=63) and cubically distributed 
patches, and 1 corresponds to a complete energy asym- 
metry (61=62=0 and e3=6tot) and a unit cell elongated in 
a single direction. The crystal symmetry group and the 
choice of A^, however, constrain 7 to within < 7 < 0.1 
and prevent dimerization [29 . 

Phase diagram - For 30 randomly selected sets of 
patch geometry and bond energies we numerically de- 
termine the solubility line using free-energy integration 
and the metastable vapor-liquid line using Gibbs Ensem- 
ble Monte Carlo simulations [29]. The gas-liquid criti- 
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FIG. 1. Temperature-density p phase diagrams of patchy models with different 7 at fixed ^ = 0-33 (left panel), and with 
different at fixed 7 = 0.0172 (right panel). Insets depicts the same phase diagrams with T rescaled following WPT. The 
crystal- fluid coexistence lines (symbols) are then close to yet distinct from the WPT solubility line (solid black line). For visual 
clarity, only the critical points (black symbols) from the gas-liquid metastable lines (symbols plus line) are reported in the 
insets. 



cal temperature Tc generally decreases with increasing 7 
because patch proximity anti-correlates bond formation 
and decreases the liquid entropy, although the limited 
number of systems studied partially hides this feature 
(Fig. [1]) [34- The solubility line, by contrast, is clearly 
similar for different geometries at fixed and monotoni- 
cally shifts to lower temperatures with increasing energy 
asymmetry (Fig. [T]) . In qualitative agreement with these 
results, the Wertheim perturbation theory (WPT) for 
associated liquids predicts that, independently of patch 
positions, /3etot + /^Mbond invariant of temperature 

along the solubility line [30l [35H37] . WPT overestimates 
the solubility temperature at all densities p, but nonethe- 
less remarkably collapses the simulation results (Fig. [l] 
insets) [29]. 

The numerical validation of WPT's Tc predictions - ac- 
curate to within 10-15% - allows us to estimate the size 
of the "crystallization gap" for a broad variety of models 
(Fig.[2| [38 . Interestingly, we find that for patch energy 
sets {ci} giving a same a lower second virial coeffi- 
cient B2 results in a larger crystallization gap (Fig. [2]A.). 
Contrary to George and Wilson's (GW) crystallization 
slot proposal that log(— B2*) < 5 identifies facile crystal- 
lization [3 , the asymmetric models reveal that B2 does 
not by itself sets the size of the crystallization gap. The 
proposal is reasonable at low but breaks down at high 

where it even includes systems for which the criti- 
cal point is fully stable (red star in Fig. |2]A. and C). In 
such a system access to the crystal from a slowly con- 
centrating, low-density solution would have to side-step 
the metastable gas-liquid coexistence regime, which typ- 
ically prevents the formation of all but the smallest crys- 
tallites [39l Ho]. High interaction asymmetry therefore 




FIG. 2. A: WPT difference between rescaled melting temper- 
ature and rescaled critical temperature T* at the critical 
density pc ~ 0.2, i.e. the "crystallization gap". Each circle 
represents a distinct {ei} realization and is colored accord- 
ing to B2{Tc). B: WPT predictions (black crosses) compared 
with simulation results for different patch geometries (sym- 
bols as in Fig. [T]). C: Phase diagram of a specific {a} for 
which WPT predicts a stable gas-liquid coexistence. Even 
for an interaction range and patch coverage that would nor- 
mally result into a metastable gas-liquid line, bond energy 
asymmetry can lift the critical point above the solubility line. 
An extrapolation of Tc using Grand Canonical Monte Carlo 
confirms this inversion [29j . 



provides a microscopic rationale for the failure of the GW 
proposal in a number of experimental systems [12 ^ • 

Crystallization - Even if crystallization is thermo- 
dynamically possible the free energy drive may be insuf- 
ficient to induce a phase transition on experimentally- 
relevant timescales. We thus examine the role of asym- 
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metry on homogeneous nucleation using umbrella sam- 
pling simulations [29 . Systems near their critical den- 
sity pc ^ 0.2 at different degrees of supersaturation 
rj = ^"^Zt ^ where Tm is the solubility temperature at 
that density, are considered. Unsurprisingly, the lower is 
the supersaturation, the higher is the free-energy barrier 
and the larger is the critical nucleus (Fig. [3^). Across 
various patch geometries qualitatively similar results are 
obtained, but increasing the energy asymmetry signifi- 
cantly lowers the chemical potential difference, /3A/i, be- 
tween the fluid and the crystal. At high bond energy 
asymmetry fewer patches dominate the energy of the 
two phases, which makes that contribution more simi- 
lar and reduces the drive to crystallize. Higher densities 
are thus needed to obtain a comparable nucleation bar- 
rier. Although this effect is not a fundamental limitation 
for particles to crystallize, proteins in high-density solu- 
tions may partially unfold and aggregate, which interferes 
with crystallization [39 . In addition, at high ( the nar- 
row crystallization gap also results in larger free-energy 
nucleation barriers. High energy asymmetry thus hinders 
nucleation kinetics. 

Classical nucleation theory (CNT) describes crystal 
reasonably well far above the critical point, but near 
and below Tc the assembly behavior is more complex. 
Previous studies of isotropically attractive systems have 
shown that well below Tc spinodal decomposition leads to 
dynamical arrest [11 , because spontaneous density fluc- 
tuations result in dense regions within which binding is 
irreversible. In similar systems near the critical point, 
"two-step" nucleation is favored [7 , because crystal for- 
mation is easier in high-density fluid regions than in low- 
density regions. The corresponding assembly behavior 
of patchy systems, whose low-density crystals may not 
be favored by spontaneous fluid density fluctuations , 
is here studied in unbiased constant NpT MC simula- 
tions. These simulations sketch out the minimum free 
energy path for the assembly, which we track along the 
largest drop and the largest crystal cluster reaction co- 
ordinates (Fig. |3]) [29]. Far above Tc, the largest cluster 
is always crystalline and CNT applies. Near the criti- 
cal point (within ~10% of Tc) a growing liquid drop first 
forms and only subsequent structural reorganization re- 
sults in a large crystal cluster. Indeed, because of the low 
nucleation barrier within a dense fluid drop (Fig. [3^) the 
drop may shelter many microcrystals. A larger crystal 
then forms from the rearrangement of these small crys- 
tals, rather than from the growth of a single critical nu- 
cleus (red versus blue symbols in Fig. [sjC), as is seen far 
above Tc. Yet, because these rearrangements are acti- 
vated, forming a large cluster can be much slower than 
classical nucleation, which is rare but rapid (blue sym- 
bols in FigjSp). Below the critical point, the bonds are 
much less labile. The system forms a large dense drop, 
but no large-scale reorganization occurs on the simula- 
tion timescale. Only small crystal clusters are thus de- 
tected. These results reveal that patchy particles have 
an assembly dynamics comparable to that of isotropic 
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FIG. 3. Rescaled size of the largest liquid cluster N* versus 
that of the largest crystal cluster N* for a model with C — ^ 
and 7 = 0.017. The liquid cluster is rescaled over the number 
of particles in the system (A^=864) and the crystal cluster 
over the size of the largest crystal cluster that fits in the sim- 
ulation box (600). The different trends correspond to starting 
homogeneous fluid configurations under different conditions, 
as illustrated in the insets on the right. Above the critical 
point, nucleation barriers can be computed as shown in inset 
B (blue circle 77 = 0.37, red squares 77 = 0.46, green trian- 
gles r] = 0.55 and magenta right-pointed triangles 77 = 0.64). 
Inset A illustrates the height of the nucleation barrier for dif- 
ferent models (blue circles C = and 7 = 0.017, red squares 
C = and 7 = 0.021, green right-pointed triangles C = 
and 7 = 0.065, black triangles C = 0.1 and 7 = 0.017, ma- 
genta left-pointed triangles ( = 0.2 and 7 = 0.017.) Inset 
C provides the size of the largest crystal cluster along the 
simulation above (blue) and at the critical point (red). 



particles, which confirms the relevance of the crystalliza- 
tion gap for identifying good crystallization conditions. 
If a large crystal is necessary for diffraction studies, one 
should target crystallization conditions well above the 
critical point, where CNT describes the dominant ho- 
mogeneous nucleation pathway [43 . Recent advances in 
X-ray laser technology, however, provide diffraction data 
from relatively small crystals [44], which should extend 
the experimentally useful range of the crystallization gap 
closer to Tc, although probably still not below Tc. 

Percolation - We finally consider direct percolation, 
which can also dynamically compete with crystalliza- 
tion. Below the percolation threshold Tp{p) the system 
forms infinitely large spanning networks that can be long- 
lasting at high bond strengths [45^. To explore the inter- 
play between percolation and bond energy asymmetry 
(patch geometry asymmetry has only a weak effect), we 
determine Tp(pc) using finite-size rescaling [29 , 46 in sys- 
tems with relatively low {( = 0.22) and high {( = 0.55) 
bond energy asymmetry (Fig. [4|. In the first system the 
percolation threshold lies just above Tc, while in the sec- 
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T MC sweeps (x 10000) 




FIG. 4. Percolation behavior for C = 0.2 (top) and C = 0.55 
(bottom). Left panels show the probability P of observing a 
spanning network as a function of temperature for system of 
size A^=2048 (blue circles), 4000 (red squares) and 6912 (green 
triangles). The superimposed vertical lines indicate the crit- 
ical temperature (dashed), the melting temperature (solid) 
and the percolation temperature (dot-dashed) estimated by 
finite-size scaling (inset). The right panels show the distribu- 
tion of bond life time in the network respectively at T=0.055 
and T=0.048 (orange stars in the left panels). Blue circles in- 
dicate the strongest bond, red squares the second strongest, 
and green triangles the weakest. 



gested m [22] . It also offers a rationalization of the GW 
crystallization slot proposal as well as for its occasional 
failure. At low bond energy asymmetry, the B2 slot pre- 
scription falls within the slot. At high asymmetry, a 
large crystallization gap is only observed for B2 below 
the slot, which corresponds to long-living gels, while for 
B2 within the slot the crystallization gap is very small or 
even negative. The GW crystallization slot is therefore a 
necessary but insufficient condition for detecting optimal 
experimental conditions. 

Although we are now closer to understanding simple, 
monomeric protein crystallization, we are still a long way 
from capturing the assembly features of more complex 
proteins. Some proteins dimerize or change conforma- 
tion on a timescale comparable to their crystallization, 
while membrane proteins typically require entirely dif- 
ferent crystallization approaches than the one considered 
in this letter. Further modifications to patchy particle 
models, such as self-interacting or dynamically evolving 
patches, may thus guide our understanding of these com- 
plex and crucial biological objects. 
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Appendix A: Model details 



ond, in which the strongest bond is much longer-lasting, 
Tp is well above the solubility line. The dynamical rele- 
vance of percolation on crystallization is estimated from 
the distribution of bond lifetimes within the crystalliza- 
tion gap. At low bond energy asymmetry, the rearrange- 
ment of all bonds is observed within a few thousand MC 
steps. At high asymmetry, by contrast, the lifetime of 
the strongest bond (blue circles) is comparable to the 
length of the simulation (10^ MC sweeps). The network 
is frozen and no rearrangement takes place. This observa- 
tion suggests that identifying the crystallization gap may 
be insufficient for crystallizing particles with high energy 
asymmetry, because a long-lasting gel caused by direct 
percolation dynamically interferes with crystallization. 

Conclusion - In order to gain insights into protein 
crystallization and soft matter assembly, we have consid- 
ered the role of patch geometry and bond energy asym- 
metry on the crystal assembly of a family of schematic 
models. We find patch geometry asymmetry to have 
a weak effect, but bond energy asymmetry to severely 
impede the crystallization thermodynamics and kinetics. 
The crystallization gap shrinks, gel formation is favored, 
and nucleation shifts to higher supersaturations. The 
union of these observations suggests that to ease iden- 
tifying the proper crystallization conditions, it is some- 
times more effective to symmetrize the directional pair 
interactions between proteins, rather than specifically 
strengthen one of them, as is sometimes implicitly sug- 



The pair- wise interaction between particles 1 and 2, 
whose centers are at distance ri2 apart, is 



(^12 



, 1^1, 1^2) = 0HS(^12) + ^[02i,2i-l(ri2, ^1, ^2) 



+02i-l,2i(^12,^l,^2)], 



(Al) 



where and Q2 are the Euler angles. A hard-sphere 
(HS) potential captures the volume exclusion 



^HS(^) 



00 r < a 
r > cr, 



(A2) 



where a is the diameter of the particle. The patch-patch 
interaction is the product of a radial and an angular com- 
ponent 

^2i,2i-l(ri2, ^1, ^2) = V^i(ri2)c^2i,2i-l(^l, ^2), (A3) 

where 



-Ci r < 
r>\i(j' 



(A4) 



^2i,2i-l(^l, ^2) 



1 0i^2i < hi and 6>2,2i-i < hi-i 







otherwise 



(A5) 
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The range of the interaction is in units of cr, is the 
semi- width of patch 2i and 0i^2i is the angle between the 
vector ri2 and the vector defining patch 2i on particle 1. 
An analogous definition holds for ^2,2i-i- For the purpose 
of this paper, and 621 are identical across patches and 
set respectively to l.la, 0.47 rad. 



Appendix B: Geometry asymmetry parameter 



The geometry asymmetry parameter is defined as 
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{h - h? + (/i - h 



2(7,2 + /|+/|) 



(Bl) 



where U represents the i^^ eigenval ue of the inertia ten- 
sor of the object represented in Fig. B.l Each patch (in 



red) carries a mass M at its center. The inertia tensor 
is computed over the collection of the weighted patches. 
The expression for 7 guarantees the its value does not 
depend on the fictitious mass M and on the distance be- 
tween the center of the patch and the center of the par- 
ticles (radius of the particle) as long as these quantities 
do not vary from one patch to the other. Patches located 
on a perfect octahedron correspond to Ii = I2 = I3, and 
consequently 7 = 0. 




of 7 ^ 0.1. Because a cubic symmetry (7 = and lim- 
iting case of P22i2i) is not realizable within the three 
screw axes symmetry of P2i2i2i, 7 is limited from below 
as well. The adopted range and width of the interaction 
and the P2i2i2i symmetry ensures that two particles 
can interact only through one bond. This condition, to- 
gether with the impossibility for a patch to interact with 
its copy, prevents dimerization. 



Appendix C: Determination of the phase diagram 

The phase diagram is determined using specialized 
Monte Carlo (MC) techniques. Gibbs Ensemble MC sim- 
ulations directly determine the coexistence densities of 
the metastable gas and liquid phases [47 . We simulate a 
total of N =1000 particles for 10^ MC cycles, where each 
cycle consists on average of N particle displacements, N 
particle rotations, N/10 particle swaps, and 2 volume 
V changes. The critical temperature Tc and density are 
then estimated using the law of rectilinear diameters |48] . 

Because the gas-liquid line is metastable, we observe 
that for certain models crystallization happens so quickly 
that we are unable to determine the gas and liquid densi- 
ties. In such cases, we estimate the critical temperature 
as the lowest temperature at which phase separation is 
not observed within 10^ MC cycles. 

To determine the fiuid-solid coexistence curve, we inte- 
grate the Clausius-Clapeyron equation starting from one 
coexistence point using a fourth-order predictor-corrector 
algorithm [5 . The coexistence point is determined using 
free energy calculations. The free energy of the fiuid is 
computed using thermodynamic integration from the free 
energy of an ideal gas [49]. The free energy of the crys- 
tal is determined using for reference the Einstein crystal 
with fixed center of mass [50] with Hamiltonian 




FIG. B.l. Sketch of the patchy particle. To determine the 
inertia tensor, we treat the particle as spherical balls (red) 
whose center is at the surface. 



The patches are randomly distributed on the surface 
under the constraint that the simply bonded crystal be- 
longs to the P2i2i2i symmetry group. As a consequence, 
patches cannot be too close to each other - otherwise 
bonded particles would overlap and the unit cell would 
stretch beyond the range of attraction A = l.la, which 
limits the achievable asymmetry and the maximum value 



N 



N 



^trans 



where f{9i,(t)i,Xi) = 1 - cos(V^i,i) + 1 - cos(V^i,2) , 
{Oi^(j)i^Xi) ^1*6 the Euler angles describing the orienta- 
tion of particle i and ifjij is the angle formed between 
the vector defining patch j of particle i and the corre- 
sponding vector in the Einstein crystal. As explained in 
Ref . [49] , the free energy of the reference Einstein crystal 
can then be written as 



^COM 
^Ein 



COM 
trans 



,COM 



where 
pa. 
and 



COM 
trans 



3N 
"2^ 



N 



In 



27V 



(CI) 



InN (C2) 



de sm{e)d(t)dx exp[-/3Sor/(6>, 0, x)] 
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TABLE B.l. Geometry parameters: the triplets of numbers represent the (x,y,z) coordinates of each patch 

patchi patch2 patchs patch4 patchs patche 7 

-0.8036 0.8036 -0.5186 -0.5186 0.3081 0.3981 

-0.5042 -0.5042 0.2731 0.2371 -0.8084 0.8084 0.0172 

-0.3163 -0.3163 0.8102 -0.8102 0.5016 0.5016 

-0.7904 0.7904 -0.4227 -0.4227 0.3571 0.3571 

-0.5184 -0.5184 0.2807 0.2897 -0.776 0.776 0.0217 

-0.3263 -0.3263 0.8617 -0.8617 0.5199 0.5911 

-0.7191 0.7191 -0.3813 -0.3813 0.3475 0.3475 

0.5659 0.5659 -0.1669 -0.1669 -0.7668 0.7668 0-0381 

-0.4032 -0.4032 -0.9093 0.9093 0.5397 0.5397 

-0.6167 0.6167 0.6103 0.6103 -0.006 -0.006 

0.6521 0.6521 -0.3099 -0.3099 -0.9579 0.9579 0-0631 

-0.441 -0.441 0.729 -0.729 0.2871 0.2871 

-0.9042 0.9042 -0.6276 -0.6276 0.2859 0.2859 

-0.3335 -0.3335 0.5386 0.5386 -0.9074 0.9074 0-0787 

-0.2669 -0.2669 0.5622 -0.5662 0.3081 0.3081 



TABLE B.2. Energy parameters 

€1 €2 €3 c 

2. 2. 2. 0.00 

1.2462 2.5482 2.2056 0.11 

2.1 2.9066 0.9934 0.21 

0.4854 2.8266 2.688 0.33 

3.5756 0.2037 2.2207 0.49 



3. 


3. 


0. 


0.50 


3.96 


1.8 


0.24 


0.55 


4.32 


1.5 


0.18 


0.64 


4.8 


0.6 


0.6 


0.79 



The calculation of atrans straightforward, but that 
of a^P^ requires either a tedious numerical integration, 
as in Ref. [51], or an analytical approximation. We opt 
for the latter using a saddle point approximation, which 



is very efficient for high values of /3Sor such as those used 
here, because the integrand is sharply peaked. Defining 
{Oo^<po^Xo) as the reference orientation in the Einstein 
crystal and changing variable a = (cos(6>), ^, x) gives 

J dO sin{0)d(j)dxQ^'p[—P'^orf{<^)] = J dcx ex.-p[— pSoj. f (ex)] ~ 

exp[-/?So,/(ao)](27r)3/2 ^ {27rf/^ 
~ (/?Sor)3/2 det(J?[/(ao)])i/2 (/?5or)3/2 det(if[/(ao)])i/2 ' 

such that 

/3a^OM « 3 j^^^^^^^ ^ 1 j^^g^ det(F[/(ao)])}, (C3) 

where det(i^[/(ao)]) is the determinant of the Hessian 
of function / computed at O-q. The analytical expression 
of /(a) is reported below defining y = cos(^) and ( as 
the angle between the vectors identifying patch 1 and 2, 
which does not depend on the orientation of the particle. 



f{oc) = \{^- 3\/l - ^2 cos(0) cos(0o) sin((9o) - \/l - ^/^ cos(x) cos(xo) sin(6>o) - 2vT^^cos(^) cos(C) sin(0o) sin(xo) sin(C) 
- 2?/cos(xo)cos(C)sin(l9o)sin(C) + 2?/cos((/))cos(x)sin((/)o)sin(x( ' ' '^'^^ ^ ' x • / x • / x • 9/>x 

-yr^ - 



^^^^v-u/ -"^v^/ : — VY^/ —v/v/ — vY^uy --^v/v^o) sin^(C) - 2?/cos((/)) cos((/)o) sin(x) sin(xo) sin^(C) 
cos(x) cos(xo) sin(^o) sin^(C) + ycos{(j)) cos(0o) cos(x) sin(^o) sin(2C) + cos(^) sin(0o) sin(x) sin(6>o) sin(2C) 
cos(^o) [—^y + cos(0o) cos(xo) sin(0) sin(x) + cos^(C) {—y + ^cos(x) cosfvn) sinfc^) sinfc^n) — cosfc^n) cosfvn) sinfc^) 



) ; tiUi[ZL^ ) cosi^^; smi^^o ; smi^x; smi^t^o ) 
cos(xo) sin(^) sin(0o) - cos(0o) cos(xo 

'v^ sin^r — 9 rnsT Vr>^ «i"n(^0Q j sin(x) sin^(^) 



os(6>o) [-3^ + cos(0o) cos(xo) sin(0) sin(x) + cos (C) {-y + ^cos(x) cos(xo) sin(^) sin(0o) - cos( 
\/l — y'^ cos(x) cos(C) sin(C) + ys\v?{Q + cos((/)o) cos(xo) sin((/)) sin(x) sin^(C) — 2 cos((/)) cos(xo! 
cos(x) cos(xo) (2cos((/)) cos((/)o) sin(C) + sin((/)) sin((/)o)(l + sin^(C))) + y'^ cos(0) cos((/)o) o 

/l - y'^ cos(xo) sin((/)) sin((/)o) sin(2C) + sin((/)) [- sin((/)o) (2 sin(x) sin(xo) sin^(C) 

^0) sin(C) + sin(xo) (-27/cos(x) sin2(C) + \/i 



\- cos 



,)sin(0) sin(x)) 
5((/)o) cos(xo) sin(2C) 
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Once the free energy of the reference crystal is known, 
the free energy of the actual crystal is obtained following 
a standard protocol [51 . Several simulations along an 
isobar starting from the fluid and from the crystal are 
then necessary to determine the temperature at which 
the chemical potential of the two phases coincides [49l 

El]. 



Appendix D: Wertheim's perturbation theory 

According to Wetheim's theory [35l[36], the free energy 
of the fluid can be approximated by the hard sphere free 
energy plus a bond free energy correction 

df = ms + ^^bond, (Dl) 



where 



/3«bond = E(l"^--^) + ?- (D2) 

aer 

Here m is the total number of attractive sites and Xa is 
the probability that the molecule is not bonded at site a. 
Similarly, the chemical potential is given by 

^IJ.f = Paf ^ /3aHS + P^^bond H \ , 

P P P 

(D3) 
(D4) 



where 



1 1 



dp J \Xa 2 



In the solid, /Sus « because the ratio ^ is 

Ps 

small [271 E2]- The energetic contribution to the free 
energy is the energy of the fully-bonded system — /3etot, 
while the entropic term is approximated using the range 
of interaction and the width of the patches [27] , 



(3i^s = PcLs = -31n(A - 1) - In 



/3etot. (D5) 



At coexistence, temperature, pressure and chemical 
potential of the fluid and the solid phase have to be 
identical. The pressure of the solid is once again ig- 
nored, so the only remaining constraint for coexistence is 
PcPf = PcPs' Using the equations above, it follows that 



PcPUS , PcPh^ 



ond 



PcdYLS + /5<^bond + 

P P 

= -31n(A - 1) - In (^^^ - /^^etot. (D6) 

As the hard- sphere system is temperature independent, 
it holds that 



<^bond + 



Pbond 



etot = C{p), (D7) 



where 

C{p) = -pam - ^ - 31n(A - 1) - In (^^^ , (D8) 

which is a constant that depends only on p. It follows 
that ttbond + ^^^^ + etot represents a good rescaling factor 
for the temperature to obtain the master solubility line 
across the different models. 



Appendix E: Grand Canonical Monte Carlo 
approximation for critical point 

To check the position of the critical point for the 
extreme case of Fig. 2C (ei = 4.6655, 62 = 1.2908 
and 63 = 0.0437), we perform Grand Canonical Monte 
Carlo (GCMC) simulations. Because one of the patches 
is markedly stronger than the others, percolation takes 
place at high temperature. Consequently, the simulation 
dynamics is slow and the estimate of the critical tem- 
perature is affected by large errors. Indeed, the GCMC 
simulations have a poor sampling. To obtain a better 
estimate of the phase diagram, we perform GCMC for 
systems with increasing strength of the strongest patch 
(keeping the other patches identical) and we fi t a p ower- 
law to the resulting critical temperatures (Fig |E.l| ). The 
value of the fit for ei = 4.6655 (T = 0.045 in units of etot) 
confirms the stability of the critical point with respect to 
the solubility line. 



^ -1.6 




0.5 1 

log(8^) 



FIG. E.l. Fit of the critical point determined with GCMC 
simulations for systems with increasing strength of the 
strongest patch. 



Appendix F: Determination of the percolation 
threshold 



To determine the temperature at which percolation 
takes places, we run 20 NVT simulations with respec- 
tively N = 2048, TV = 4000 and N = 6912 at several 
temperatures and p = 0.2. During the simulation, we 
determine the size of the biggest network defined as the 
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largest set of particles connected by at least two bonds. If 
such a cluster grows sufficiently to span the whole simula- 
tion box along one dimension within 10^ MC sweeps, the 
system is deemed percolating. The probability of perco- 
lation is the fraction of simulations that show a percolat- 
ing cluster. For finite-size rescaling we use the tabulated 
critical exponents and the standard procedure [46 . The 
data are in agreement with the tabulated 3D critical ex- 
ponent to within 1%. 

Appendix G: Definition of crystal cluster and 
umbrella sampling simulations 

We determine the size of the crystal clusters in the sim- 
ulation box following a standard procedure that defines 
crystal-like bond and crystal-like particles [53]. Due to 
the highly specific patch-patch interaction of our model, 
we straightforwardly define a crystal-like bond between 
particles 1 and 2 when they are bonded in our model: 
r < Act, 6i^2i < ^ and 6>2,2i-i < ^ for some i. For 
C > 0.33, the weakest interactions tend to be looser, so we 



relax the conditions and define a bond when r < 1.2<j, 
^i,2i < 2(5 and ^2,2i-i < 2^. A particle is considered 
to be crystal-like if it has six crystal-like bonds. Two 
crystal-like particles belong to the same crystal cluster if 
a crystal-like bond connects them. Visual inspection of 
these "crystals" confirms that the criterion selects actual 
crystal clusters. 

For the umbrella sampling simulations, we use a bias- 
ing harmonic potential with spring constant that varies 
between 0.07 and 0.12, depending on the model and the 
temperature 

H^'^' = f^{n-no)^, (Gl) 

where n represents the size of the largest crystal clus- 
ter and no is the cluster size corresponding to sampling 
window. The sampling windows are typically positioned 
every 3 particles, although sometimes denser sampling 
is required. The results of each simulation are then an- 
alyzed following the umbrella sampling standard proto- 
col [53j . 
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